clear all; clc; close

load EcoliHyperShockEpi.mat
i = 6;
     figure
     hold on
     
     xlabel('Time (s)')
     ylabel('A/A0')
     title(experiment{i})
     
     
     [numCells,~]=size(cell_width{i}); % numCells is number of cells
%      R = 0.000555556; % The average unperturbed growth rate
%      deltaT = R.*(t_vec{i});
%      L0 = exp(deltaT);
   
     [~,numPoints]= size(t_vec{i});
     
     locationOfShock = 32;
     
     %numPoints = locationOfShock; % when plotting only before the jump. 
 
     
    %% calculating the average width as function of time
     avg_width = zeros(numPoints,1);
     N = zeros(numPoints,1);
     for j=1:numCells
         if (~isnan(cell_width{i}(j,locationOfShock))) && (~(cell_width{i}(j,locationOfShock)==0))
           
           
           plot(t_vec{i},cell_width{i}(j, :)/cell_width{i}(j, locationOfShock)) %/cell_width{i}(j, locationOfShock)
           
           normalized_width = cell_width{i}(j,:)/cell_width{i}(j,locationOfShock);
           
           normalized_width(isnan(normalized_width)) = 10;

        for k=1:numPoints
            if normalized_width(k) ~= 10
           
         avg_width(k) = (normalized_width(k) +  avg_width(k));
         N(k) = N(k) + 1;
            end
        end
         end
       % plot(deltaT,normalized_width)
     end
    
     avg_width = avg_width ./ N;
   
     
     figure
     hold on
     
     plot(t_vec{i}, avg_width)
     xlabel('Time (s)')
     ylabel('A/A0')
     title(experiment{i})
     
     fig2pretty
    

     
        %% calculating the standard deviation corresponding to the mean above
     std_width = zeros(numPoints,1);
     N = zeros(numPoints,1);
     for j=1:numCells
         if (~isnan(cell_width{i}(j,locationOfShock))) && (~(cell_width{i}(j,locationOfShock)==0))
           
           normalized_width = cell_width{i}(j,:)/cell_width{i}(j,locationOfShock);
         
           normalized_width(isnan(normalized_width)) = 10;
       
        for k=1:numPoints
            if normalized_width(k) ~= 10
           
         std_width(k) = std_width(k) +  (normalized_width(k) -  avg_width(k))^2;
         N(k) = N(k) + 1;
            end
        end
         end
      
     end

     std_width = sqrt(std_width ./ N);

     
%      figure
%      hold on
%      
%      plot(t_vec{i}, std_width)
%      xlabel('Time (s)')
%      ylabel('std')
%      title(experiment{i})
%      
%      fig2pretty
     
     
     
%      figure
%      hold on
%      
%      plot(t_vec{i}, avg_width, 'b-')
%      
%      plot(t_vec{i}, avg_width - 2*std_width, 'r-')
%      
%      plot(t_vec{i}, avg_width + 2*std_width, 'r-')
%      
%      xlabel('Time (s)')
%      ylabel('A/A0')
%      title(experiment{i})
%      
%      fig2pretty
    
     
     %% Export arrays
       
     experiment2 = experiment;
     
     experiment2{1}(3:4) = 'to';
     experiment2{2}(4:5) =  'to';
     experiment2{3}(3:4) =  'to';
     experiment2{4}(3:4) = 'to';
     experiment2{5}(3:4) =  'to';
     experiment2{6}(3:4) =  'to';
     
     writematrix(t_vec{i},['time_'  experiment2{i} '.csv']) 
     writematrix(avg_width,['avgwidth_' experiment2{i} '.csv']) 
     writematrix(std_width,['stdwidth_' experiment2{i} '.csv']) 
 